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Abstract 

A non-overlapping domain decomposition algorithm is proposed to 
solve the linear system arising from mixed finite element approximation 
of incompressible Stokes equations. A continuous finite element space for 
the pressure is used. In the proposed algorithm, Lagrange multipliers 
are used to enforce continuity of the velocity component across the sub- 
domain domain boundary. The continuity of the pressure component is 
enforced in the primal form, i.e., neighboring subdomains share the same 
pressure degrees of freedom on the subdomain interface and no Lagrange 
multipliers are needed. After eliminating all velocity variables and the in- 
dependent subdomain interior parts of the pressures, a symmetric positive 
semi-definite linear system for the subdomain boundary pressures and the 
Lagrange multipliers is formed and solved by a preconditioned conjugate 
gradient method. A lumped preconditioner is studied and the condition 
number bound of the preconditioned operator is proved to be independent 
of the number of subdomains for fixed subdomain problem size. Numerical 
experiments demonstrate the convergence rate of the proposed algorithm. 

keywords domain decomposition, incompressible Stokes, FETI-DP, BDDC 
AMS 65F10, 65N30, 65N55 

1 Introduction 



Domain decomposition methods have been studied well for solving incompress- 
ible Stokes equations and similar saddle-point problems; see, e.g., [TBI 1231 1201 
[T01 [31 [221 [HI 12H1 HH 125] • In many of those work, special care need be taken to 
deal with the divergence-free constraints across subdomain boundaries, which 
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often lead to large coarse level problems. The large coarse level problem will be 
a bottleneck in large scale parallel computations, and additional efforts in the 
algorithm are needed to reduce its impact, cf. [33 1321 1301 El 111 HS1 133] ■ Some 
recent progress has been made by Dohrmann and Widlund [H |6] for the almost 
incompressible elasticity, where the coarse level space is built from discrete sub- 
domain saddle-point harmonic extensions of certain subdomain interface cut-off 
functions and its dimension is much smaller than those in the previous studies. 
Kim and Lee [El H3 H2 with Park] studied both the FETI-DP and BDDC 
algorithms for incompressible Stokes equations where a lumped preconditioner 
is used and reduction in the dimension of the coarse level space is also achieved. 

In most above mentioned applications and analysis of domain decompo- 
sition methods for incompressible Stokes equations, the mixed finite element 
space contains discontinuous pressures. Application of discontinuous pressures 
in domain decomposition methods is natural. The decomposing of the pressure 
components to independent subdomains can be handled conveniently and no 
continuity of pressures across the subdomain boundary need be enforced. How- 
ever, a big class of mixed finite elements used for solving incompressible Stokes 
and Navier-Stokes equations have continuous pressures, e.g., the well known 
Taylor- Hood type [27]. There have been a variety of approaches using con- 
tinuous pressures in domain decomposition methods for solving incompressible 
Stokes equations, e.g., by Goldfeld [9], by Sistek et. al. [26], and by Benhassine 
and Bendali 1 . In their work, an indefinite system of linear equations need be 
solved, either by a generalized minimal residual method or simply by a conju- 
gate gradient method. To the best of our knowledge, no scalable convergence 
rate has been proved analytically for any of those approaches using continuous 
pressures. 

In this paper, we propose a non-overlapping domain decomposition algo- 
rithm for solving incompressible Stokes equations with continuous pressure fi- 
nite element space. The scalability of its convergence rate is proved. In this 
algorithm, the subdomain boundary velocities are dealt with in the same way 
as in the FETI-DP method: a few for each subdomain are selected as the coarse 
level primal variables, which are shared by neighboring subdomains; the others 
are subdomain independent and Lagrange multipliers are used to enforce their 
continuity. The subdomain boundary pressure degrees of freedom are all in the 
primal form. They are shared by neighboring subdomains and no Lagrange mul- 
tipliers are needed for their continuity. After eliminating all velocity variables 
and the independent subdomain interior parts of the pressures, the system for 
the subdomain boundary pressures and the Lagrange multipliers is shown to be 
symmetric positive semi-definite. A preconditioned conjugate gradient method 
with a lumped preconditioner is studied. As strong condition number bounds 
as for the scalar elliptic case are established. In the proposed algorithm and 
in the estimate of its condition number bound, no additional coarse level vari- 
ables, except those necessary for solving scalar elliptic problems, are required 
for incompressible Stokes problems. The resulting coarse level problem is also 
symmetric positive definite. 

To stay focused on the purpose of this paper, the discussion of the proposed 
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algorithm and its analysis are based on two-dimensional problems, even though 
the same approach can be extended to the three-dimensional case without sub- 
stantial obstacles. It is also worth pointing out that the domain decomposition 
algorithm and its analysis presented in this paper apply equally well, with only 
minor modifications, to the case where discontinuous pressures are used in the 
mixed finite element space. 

The remainder of this paper is organized as follows. The finite element dis- 
cretization of the incompressible Stokes equation is introduced in Section O A 
domain decomposition approach is described in Section |3l The system for the 
subdomain boundary pressures and the Lagrange multipliers is derived in Sec- 
tion 2J Section [5] provides some techniques used in the condition number bound 
estimate. In Section [6j a lumped preconditioner is proposed and a scalable 
condition number bound of the preconditioned operator is established. At the 
end, in Section [3 numerical results for solving a two-dimensional incompressible 
Stokes problem are shown to demonstrate the convergence rate of the proposed 
algorithm. 



2 Finite element discretization 

We consider solving the following incompressible Stokes problem on a bounded, 
two-dimensional polygonal domain with a Dirichlet boundary condition, 

!-Au + Vp = f , in tt , 
-V • u = 0, in il , 
u = u on , on dQ , 

where the boundary data 11,90 satisfies the compatibility condition J gn u^o ■ n = 
0. For simplicity, we assume that uqq = without losing any generality. 

The weak solution of {TJ is given by: find u £ (i?^(0)) 2 = {v £ {H x {9)) 2 | v = 
on dil} and p £ L 2 (il), such that 




a(u,v) + 6(v,p) = (f,v), Vve (ffoHfi)) 2 , 
b(u,q) = 0, Vgei 2 (0), 



where 

a(u,v)= f Vu-Vv, b(u,q) = - /(V-u)g, (f,v)=/fv. 

JO JO Jo 

We note that the solution of ([2]) is not unique, with the pressure p different up 
to an additive constant. 

A modified Taylor-Hood mixed finite element is used in this paper to solve 
(0 ■ The domain f2 is triangulated into shape- regular elements of characteristic 
size h. The pressure finite element space, Q C L 2 (il), is taken as the space of 
continuous piecewise linear functions on the triangulation. The velocity finite 
element space, W £ (Hq(Q)) , is formed by the continuous piecewise linear 
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functions on the finer triangulation obtained by dividing each triangle into four 
subtriangles by connecting the middle points of its edges. A demonstration of 
this mixed finite element on a triangulation of a square domain is shown in 
Figure [TJ 



Figure 1: A modified Taylor- Hood mixed finite element 
The finite element solution (u,p) G W0Q of ([2|) satisfies 



(3) 
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where A, B, and f represent respectively the restrictions of a(-, •), &(•, •) and 
(f , ■) to the finite-dimensional spaces W and Q. We use the same notation in 
this paper to represent both a finite element function and the vector of its nodal 
values. 

The coefficient matrix in ([3]) is rank deficient. A is symmetric positive def- 
inite. The kernel of B T , denoted by Ker(B T ), is the space of all constant 
pressures in Q. The range of B, denoted by Im(B), is orthogonal to Ker(B T ) 
and is the subspace of Q consisting of all vectors with zero average. The solution 
of ([3]) always exists and is uniquely determined when the pressure is considered 
in the quotient space Q/Ker(B T ). In this paper, when q 6 Q/Ker(B T ), q 
always has zero average. For a more general right-hand side vector (f , g) given 
in ([3]), the existence of its solution requires that g S Im(B) 7 i.e., g has zero 
average. 

The modified Taylor-Hood mixed finite element space W x Q, as shown in 
Figure [H is inf-sup stable in the sense that there exists a positive constant /3, 
independent of h, such that 




D velocity and pressure 



• velocity only 



(4) 




>/3|MIl2, Vq e Q/ Ker(B T ), 



cf. [21 Chapter III, §7], or equivalently in matrix/vector form, 



(5) 



W 



( 5 ,Bw) 2 

iup 1 FT 

eW (w, Aw) 



>f3 2 (q,Zq), VqeQ/Ker(B T ). 
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Here, as always in this paper, (•, •) represents the inner product of two vectors. 
The matrix Z represents the mass matrix defined on the pressure finite element 
space Q, i.e., for any q G Q, \\q\\ 2 L2 = (q, Zq). It is easy to see, cf. jMJ Lemma 
B.31], that Z is spectrally equivalent to h 2 I for two-dimensional problems, where 
I represents the identity matrix of the same dimension, i.e., there exist positive 
constants c and C, such that 

(6) ch 2 I <Z< Ch 2 I. 

Here, as in other places of this paper, c and C represent generic positive con- 
stants which are independent of the mesh size h and the subdomain diameter 
H (discussed in the following section). 



3 A non-overlapping domain decomposition ap- 
proach 

The domain fl is decomposed into N non-overlapping polygonal subdomains flj , 
i = 1,2, N. Each subdomain is the union of a bounded number of elements, 
with the diameter of the subdomain in the order of H. The nodes on the 
interface of neighboring subdomains match across the subdomain boundaries 
T = (U<9f2i)\9f2. T is composed of subdomain edges, which are regarded as 
open subsets of T, and of the subdomain vertices, which are end points of edges. 

The velocity and pressure finite element spaces W and Q are decomposed 
into 

W = W,0W r , Q = Q/0Qr, 

where W/ and Qi are direct sums of independent subdomain interior velocity 
spaces Wj\ and interior pressure spaces Qj, respectively, i.e., 

N N 

w, = 0wf, Qi = 0Qf. 

i=l i=l 

Wr and Qr are subdomain boundary velocity and pressure spaces, respectively. 
All functions in Wr and Qr are continuous across the subdomain boundaries 
T; their degrees of freedom are shared by neighboring subdomains. 

To formulate our domain decomposition algorithm, we introduce a partially 
sub-assembled subdomain boundary velocity space Wr, 



W r = W n 0W A = Wn0 0W; 




Here, Wn is the continuous, coarse level, primal velocity space which is typi- 
cally spanned by subdomain vertex nodal basis functions, and/or by interface 
edge basis functions with constant values, or with values of positive weights on 
these edges. The primal, coarse level velocity degrees of freedom are shared by 
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neighboring subdomains. The complimentary space Wa is the direct sum of in- 

(i) 

dependent subdomain dual interface velocity spaces W A , which correspond to 
the remaining subdomain boundary velocity degrees of freedom and are spanned 
by basis functions which vanish at the primal degrees of freedom. Thus, an el- 
ement in the space Wr typically has a continuous primal velocity component 
and a discontinuous dual velocity component. 

The functions wa in Wa are in general not continuous across T. To enforce 
their continuity, we define a boolean matrix Ba constructed from {0, 1, —1}. On 
each row of -Ba, there are only two non-zero entries, 1 and —1, corresponding 
to the same velocity degree of freedom on each subdomain boundary node, 
but attributed to two neighboring subdomains, such that for any wa in Wa, 
each row of -Ba w a = implies that these two degrees of freedom from the 
two neighboring subdomains be the same. When non-redundant continuity 
constraints are enforced, Ba has full row rank. We denote the range of Ba 
applied on Wa by A, the vector space of the Lagrange multipliers. 

In order to define a certain subdomain boundary scaling operator, we intro- 
duce a positive scaling factor 5>(x) for each node x on the subdomain bound- 
ary r. Let N x be the number of subdomains sharing x, and we simply take 
S'(x) = 1/N X . In applications, these scaling factors will depend on the heat 
conduction coefficient and the first of the Lame parameters for scalar elliptic 
problems and the equations of linear elasticity, respectively; see [TOUTS] , Given 
such scaling factors at the subdomain boundary nodes, we can define a scaled 
operator Ba,d- We note that each row of Ba has only two nonzero entries, 1 
and —1, corresponding to the same subdomain boundary node x. Multiplying 
each entry by the scaling factor S^(x) gives us Ba : d- 

Solving the original fully assembled linear system ([3]) is then equivalent to: 
Pr, A) e W ; © © W A W n g r © A, such that 



find(u/, p/, u A , u n 

A n 
B u 
Aai 
Am 

Bn 




(7) 
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where the sub-blocks in the coefficient matrix represent the restrictions of A 
and B in to appropriate subspaces. The leading three-by-three block can 
be made block diagonal with each diagonal block representing one independent 
subdomain problem. 

Corresponding to the one-dimensional null space of ([3]) , we consider a vector 
of the form (u/, p/, ua, u n , pr, A) = (0, 1 PI , 0, 0, l Pr ,A), where l Pl e Qi 
and l Pr £ Qr represent vectors with value 1 on each entry. Substituting it into 
([7]) gives zero blocks on the right-hand side, except at the third block 



(8) 



fA - [BJ a B ta] 



L P1 



B T A \. 
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The first term on the right-hand side represents the line integral of the normal 
component of the velocity finite element basis functions across the subdomain 
boundary on neighboring subdomains. Corresponding to the same subdomain 
boundary velocity degree of freedom, their values on the two neighboring sub- 
domains are negative of each other. Therefore 



in 



Pr 



[bJ a B- A ] 

from which we know that fA = 0, for 

A = —Ba,d[BJ a B- a ] 



b a b a , d [bJ a b- a ] 



l>T 



IT 



Therefore, a basis of the one-dimensional null space of is 



(9) 



0, l Pl , 0, 0, l Pr , -B A , D [Bf A Bf A ] 



4 A reduced symmetric positive semi-definite 
system 

The system ([7]) can be reduced to a Schur complement problem for the variables 
{pr, A). Since the leading four-by-four block of the coefficient matrix in (0 is 
invertible, the variables (uj, pj, u Al un) can be eliminated and we obtain 



(10) 

where 
(11) 

G 



G 
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(12) 9 
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We denote 



(13) 



A = 



An 




AlA 


A m 


Bji 





B IA 


Bm 


A A i 


b Ja 


A A a 


A A n 


Am 


B m 


A nA 


Ann 



and Bn = 



Bn B F a B T n 
Ba 



We can see that — G is the Schur complement of the coefficient matrix of ([7]l 
with respect to the last two row blocks, i.e., 



I 



A 

B c 



B T c 




-A-'Bl 
I 





'a 




-G 



From the Sylvester's law of inertia, namely, the number of positive, negative, 
and zero eigenvalues of a symmetric matrix is invariant under a change of co- 
ordinates, we can see that the number of zero eigenvalues of G is the same as 
the number of zero eigenvalues (with multiplicity counted) of the original coef- 
ficient matrix of (O, which is one, and all other eigenvalues of G are positive. 
Therefore G is symmetric positive semi-definite. The null space of G is derived 
from the null space of the original coefficient matrix of ([7]), and its basis is given 
by, cf. ®, 



l P r: -Ba,d[BJa b y a] 



1 



PT 



We denote X = Qr©A. The range of G, denoted by Rq, is the subspace 
of X orthogonal to the null space of G, and has the form 



(14) R G = 



9 P v 
fix 



G X 
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The restriction of G to its range Rg is positive definite. The fact that 
the solution of ([7]) always exists for any given (f/, fA, fn) on the right-hand 
side implies that the solution of (fTo|) exits for any g defined by (fT2|). Therefore 
g G Rg- When the conjugate gradient method (CG) is applied to solve (|10l) with 
zero initial guess, all the iterates are in the Krylov subspace generated by G and 
<7, which is also a subspace of Rg, and where the CG cannot break down. After 
obtaining (pr, A) from solving (fTQ|) . the other components (uj, pi, ua, un) in 
([7| are obtained by back substitution. 

In the rest of this section, we discuss the implementation of multiplying G 
by a vector. The main operation is the product of A~ x with a vector, cf. (jTTJ) 
and (HU). We denote 



A = 



A n Bfj A IA 
Bn Bja 
Aai Bf A Aaa 



, An r — A^n — [Am Bj n AnA] , fr — 



f/ 



fA 
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and define the Schur complement 

Sn = Ann — Au r A r ^A r u, 

which is symmetric positive definite from the Sylvester's law of inertia. Sri 
defines the coarse level problem in the algorithm. The product 

" A n Bjj A IA A m " 

Bji Bia Bm 

A A i Bj A A aa A An 

. Am Bj n An A Ann . 

can then be represented by 

•^n 1 (fti - An r A r rf r ) , 

which requires solving the coarse level problem once and independent subdomain 
Stokes problems with Neumann type boundary conditions twice. 





fA 

. fn . 



A- 1 f 




+ 



rr 



5 Some techniques 

We first define certain norms for several vector /function spaces. We denote 
(15) W = W/0W r . 

For any w in W, we denote its restriction to subdomain fij by w''l A 
subdomain-wise iJ 1 -seminorm can be defined for functions in W by 

N 

|w|^=]T|w«|| 1(ni) . 
i=l 

We also define 

w = w J 0g / 0w A 0w n , 

and its subspace 
(16 i 

W = I w = (wj, pi, wa, wn) e W | B/jw/ + Bi A w A + B ;n w n = [ . 
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For any w — (w/, pi, wa, wn) G Wo, let w = (w/, wa, wn) G W. Then 
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i.e., (-, ■) 2 defines an inner product on Wq. In (|17p . the superscript W is used to 
represent the restrictions of corresponding vectors and matrices to subdomain 

Since W is essentially the subspace of W with continuous subdomain bound- 
ary velocities, the inf-sup condition ^ and (0 also holds for the mixed space 
W x Q. Denote 



(18) B 
as in (O, then 
(19) 



Bn Bja Bin 
Bti Bya -Bra 



.4 



q, Bw 

sup ' = \ > (3 2 (q, Zq) , 
wew ( w, Aw 



An A IA A m 
Aai Aaa AAn 
Am AnA Ann 



Vq G Q/Ker{B T ), 



where (3 is the same as in (0} and J3]). 

We also have the following lemma on the stability of the operator B. 



Lemma 1 For any w G W and q G Q, (Bw, q) < |w|#i ||<z||l 



Proof: 




The finite element space for subdomain boundary pressures, Qr, is a sub- 
space of L 2 (T). For each pr G Qr, its finite element extension by zero to the 
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interior of subdomains is denoted by p® , which equals pr on all subdomain 
boundary nodes and equals zero on all subdomain interior nodes. We can see 
that p^ £ Q C L 2 (Q), and \\p^ \\ 2 L 2^ = (p®,p®) z , from the definition of Z in 
Section H 

From (|11[) and (fT5j). we can see that 

G = BcA~ 1 Bq. 

In particular, we denote the first row of Bq by 

Br = [Br i Br a Brn] 

for the second row, we denote the restriction operator from W onto Wa by Ra , 
such that for any w = (wj, pi, wa, wn) £ W, Raw — wa- Then G can be 
represented by the following two-by-two block structure 



(20) G = 

where 



Gp r p r Gp r \ 
G\ Pr G\\ 



Cprpr — Br A 1 i?p , G Pr \ — BrA 1 R\B^ : 
G\ Pr = BaRaA^B^ , G\\ = BaRaA~ 1 R a ^B^. 

The pressure components of all vectors in Rq with = 0, cf. (1141) . form 
a subspace of Qr and we denote this subspace by Rq\q t . From the definition 
of we can see that for any vector pr £ Rc\Q r > Pr^vr — 0j an d then its 
extension by zero to the interior of subdomains, p^ , also has zero average. 

The following lemma follows essentially from [531 Lemma 9.1]. 

Lemma 2 For all pr £ Rg\q f , 

P 2 \\Pr\\h m < (pr,Gp r p r pr) < \\p^\\h m , 

where p^ represents the extension by zero of pr to the interior of subdomains, 
and j3 is the same as in and f5|). 

Proof: Note that even though^! -1 is indefinite in W, it is positive definite 
when restricted to a subspace of W, where the pressure component equals zero, 
and the norm || ■ is well defined. 

To prove the left side inequality, denote for any v = (vj, va, vn) £ W, 
v 1 * = (v/, 0, va, vn) £ W. We have 

(v^B^pX (p^BrA- 1 ^ 
pr,B T A- l BTp T ) = [[B^prf^ = sup X 2 ' A ~ x = sup V +r ~ - - 

SU P ' iTj = SU P > P 2 (Pr,Pr) z = P 2 \\Pr\\l H n), 

wew wT A ™ J wew w T 4w 
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where we have used the inf-sup condition (TT§)) for the inequality in the middle. 

To prove the right side inequality, for any given pr 6 Ra\Q r , denote = 
(vr, pi, va, vn) = A^ 1 B- pr, and the shorter vector v = (vj, va, Vn)- From 
the continuity of B in Lemma [T] and fj 1 7[) . we have 

(jjr.iMr^rW) = (pr.iW) = (p^,Bv) < \\p$\\ L * |v| ff i 

= |bf?|U 2 J (A-iBg-pr, I-i^pr) . = HpFII^ (pt, BrA^B^r)^ . □ 

The following corollary of Lemma [2] is an immediate result from ([6]) and the 
facts that ||pF||| 2(n) = <pF,pF) z , (pF,Pr) = (pr,Pr)- 

Corollary 1 There exist positive constants c and C , such that 

ch 2 ft 2 I Pr < G PrPr < Ch?I Pv 

where I Pv is the identity matrix of the same dimension as G PrPr , and ft is the 
same as in and 

Remark 1 Lemma\^and Corollary\]\ are not used in our proof of the condition 
number bound in Section [SI However, it is intuitive to see from Corollary [7] 
that the first diagonal block G PrPr in matrix G can be approximated spectrally 
equivalently by the identity matrix multiplied by h 2 , which is what is being done 
in our block diagonal preconditioner discussed in Section^ 

We also need define ja^ certain jump operator across the subdomain bound- 
aries T. Let P D : W -> W, be defined by, cf. [21], 

Pd = R A B AD B A R A . 

We can see that application of Pd to a vector essentially computes the difference 
(jump) of the dual velocity component across the subdomain boundaries and 
then distributes the jump to neighboring subdomains according to the scaling 
factor 6'(x). In fact, the dual velocity component is the only part of the vector 
involved in the application of Pd', all other components are kept zero and are 
added into the definition to make Pd more convenient to use in the presentation 
of the algorithm. We also have, for any w = (w_r, pi, w A , Wpi) S W, 

(P D w,P D w)^ = (B A D B A w A , B A D B A w A ) . 
The following lemma can be found essentially from [23, Section 6]; see also ([T7|) . 

Lemma 3 There exists a function Q(H/h), such that for all w € Wq, 
(P d w, P d w) x < ®{H/h) (w, w) A . 
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Remark 2 Just as for the positive definite elliptic problems discussed in ]2Si 
Section 6], for two-dimensional problems, when only subdomain corner velocities 
are chosen as coarse level primal variables, $(H/h) — C(H/h) (1 + log (H/h)); 
when both subdomain corner and edge-average velocity degrees of freedom are 
chosen as primal variables, <f>(H/h) = CH/h. 



The following lemma is also used and can be found at |10| Lemma 2.3]. 

Lemma 4 Consider the saddle point problem: find (u,p) G Q, such that 



(21) 



" A 


B T ' 




u 




' f " 


B 







P 




_ 9 _ 



where A and B are as in {3]), f G W. and g G Im(B) C Q. Let f3 be the inf-sup 
constant specified in {5|). Then 



where Z is the mass matrix defined in Section^ 

6 A lumped preconditioner 

The lumped preconditioner was first used in the FETI algorithm [7] for solving 
positive definite elliptic problems. Compared with the Dirichlet preconditioner, 
also used for the FETI algorithm [S], the lumped preconditioner is less effec- 
tive in the improvement of convergence rate, but it is also less expensive in the 
computational costs. The main operation in the lumped preconditioner is sub- 
domain matrix and vector products, while the implementation of the Dirichlet 
preconditioner requires solving subdomain systems of equations. In this pa- 
per, we discuss only the lumped preconditioner in our algorithm for solving the 
incompressible Stokes equation; study of the Dirichlet preconditioner will be 
addressed in forthcoming work. 

We consider a block diagonal preconditioner for (TTU1) . From Corollary [TJ the 
inverse of the first diagonal block G PrPr of G can be effectively approximated 
by 1/h 2 times the identity matrix. The inverse of the second diagonal block 
can be approximated by the following lumped block 




l|u|U< ||f|U-i+ flIMIz-i, 



M~ 



-l 



Ba,dRaARaB^ 



This leads to the lumped preconditioner 



M 



-l 



-l 



for solving (|10l) . 
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Remark 3 The mesh size h is used in the above preconditioner. For applica- 
tions where the mesh size is not explicitly provided and only the coefficient matrix 
in is given, an estimate of h can be obtained by comparing the nonzero en- 
tries in A and B blocks. From the definition of A and B for the incompressible 
Stokes problem @) ; entries in A and entries in B have a difference of factor h 
in general. 

M _1 is symmetric positive definite. Multiplication of M" 1 by a vector 
requires mainly the product of A with a vector. When the CG iteration is 
applied to solve the preconditioned system 



(22) M~ l G 



Pr 
A 



= M- l 9l 



with zero initial guess, all the iterates belong to the Krylov subspace generated 
by the operator M~ 1 G and the vector M~ 1 g 1 which is also a subspace of the 
range of M~ 1 G. We denote the range of M~ 1 G by Rm~ 1 g- The following 
lemma shows that the CG iteration applied to solving (|22|) cannot break down. 



Lemma 5 Let the preconditioner M 1 be symmetric positive definite. The CG 
iteration applied to solving with zero initial guess cannot break down. 



Proof: We just need to show that for any ^ x E Rm- 1 Gi Gx ^ 0. Let 
/ i = M~ 1 Gy, for a certain y 6 X and y ^ 0. Gx = GM~ 1 Gy, which cannot 
be zero since Gy ^ and y T GM~ 1 Gy ^ 0. □ 

Lemma 6 Let M^ 1 be symmetric positive definite. For any x — (pr, A) G 
Rm- 1 Gj 

(Mx,x) = max 



yeRcy^o (M^y^y) ' 
Proof: Denote the range of M~2G by R M - i/2 G . For any x s Rm- 1 Gi 

2 



(Mx,x) = (m%x,M%x\ 



max 



M?x,z 



*eR M - 1/2a ,z^o (z,z) 
(M^x,M-iy) 2 _ {y , x ) 



= max r- = max 7— — : . □ 

yeRcy^o / M~iy, M'^y) v^o-y^ {M^y, y) 



In the following, we establish a condition number bound of the precondi- 
tioned operator M _1 G. We first have the following lemma. 

Lemma 7 For any w € Wq, 

(M^Bcw, Bcw) < <&(H/h) IAw,w 



where <fr(H/h) is as defined in Lemma\^ 
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Proof: Given w = (w/, q I} wa- w n ) £ Wo, let g Pr = B r /W/ + 5 F awa 
Bmwn. We have 



(M- l B G w,Bcw) 



1 / ~ \ T ~ 



h 2 



72 (9pv9pt) + (BaBaw) Ba.dRaAPl^B^ d (^BaRawJ 



(23) 



/; 2 
1 



(9pri9pr) + (PdW,P D w) A 



{9pv,9 P v) + $(H/h) (w,w) x , 



where we used Lemma[3]for the last inequality. It is sufficient to bound the first 
term of the right-hand side in the above inequality. 

We denote w = (w;, wa, wn) £ W. Since B//w/ + B/awa + B/nwn = 0, 
cf. we have 



= ( Bw, Bw 



B//w/ + B/awa + B/nWn 
Bnwi + BrAWA + -Brawn 



-B//w/ + %wa + -B/nwn 
Br/w/ + BrAWA + Brawn 



where B is defined in (IT8l) . From ([6]) and the stability of B, cf. Lemma [T] we 
have 



1 



|w|gl||g|||j ~, |2 ^/ \ 

< Cmax jf-TTo — — = C\w\ H i = C (w,w)7, 

qeQ \\q\\ L 2 

where for the last equality, we used the fact that B//W/ +B/awa +B/nwn = 0, 
and (HZJ). □ 

Lemma 8 For any given y = (g Pr ,g\) £ Be there exits w £ Wo, such that 
B c w = y, and ^Aw,w^ < jj? (M~ 1 y,y). 

Proof: Given y = (g Pr ,gx) £ Rg, take w^ = B^ D g x . Let w^ = 

(0, w^OjeWjeWA^WnandwW = (0, 0, w£ } , 0) £ W/ Q/ W A W n . 

We have 



Bw, q 



(25) 
and 

(26) B c w [I) 



wW|^=<A AA w^,w^>, 



a) „ r co 



Br/ BrA Brn 
B A 







-Ba,d9a 




BrAW^ 
,9a 
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where we used the fact that BaB^ d = I. 

We consider the solution to the following fully assembled system of linear 
equations of the form©: find (wj"\ q ( j U \ w["\ q^ H) ) £ Wj Qi W r Qr, 
such that 



(27) 



" A H 




Air 


B TI 


B H 





Bit 





An 


Bj r 


Att 




. Bn 





Btt 






w} ' 




' 


(II) 
Qi 













(II) 







where a particular right-hand side is chosen. We first note that, since (g Pr , g\) £ 
Rc, the right-hand side vector of the above system satisfies, cf. ((HJ), 



(-B/Awi /) ) T l w + ( ffpr -SrAW^) T l Pr = gTl Pr -gl £ a ,d (Sf A l Pi + Sf A l 



i.e., it has zero average, which implies existence of the solution to (|27|) . 

Denote w^ 77 ) = (wj , w^' 1 ) <E W. From the inf-sup stability of the 
original problem ([3]) and Lemma HI we have 
(28) 

,J~, -I 2 r- ,T\ -I 2 





1 



-BjAwi J) 



< 



1 



_ BrAW^ 



1 

'^2 



5pr 



The first term on the right-hand side of 
way as done in (I24[) , and we have 



can be bounded in the same 



(29) 



S/aw^ 
BrAW^' 



<C(A AA wi /) ,wi /) 



z- 1 



the second term can be bounded by, using ([6]) 

2 

i! 

(30) 





fJlT 



c 



z- 1 



— ^2 \9pr ' $Pr ) 



Split the continuous subdomain boundary velocity w{! 7 ' into the dual part 

w A S Wa and the primal part w^ 7 ^ £ Wn, and denote i// 77 ' = (w/ , q\ IT \ w^'', 
We have, from (l27l). 



(31) 



[ fl/i B/a -B/n ] 



r (u) 



1i 



(ii) 



in) 

r (U) 



-B/Awi r) , 
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and 



(32) B c w [II) = 



Br i BrA Brn 
B A 



(n) 
w} ' 

(ii) 
1i 

(ii) 

W A 

„,(H) 



g Pr - BrA^A 




Let w = + w( n \ We can see from (|3"Tj) that w G Wq, cf. ([T 
also see from (|26| and (|32|) that Bcw = y. Furthermore, by (IT7|) . 



C 



Ml = |w«+w<">|k < |w«| 2 ffl + |w(")| 2 H1 < ^ (A AA w«,w« 

where we used (|25|), (|28]l. (|29|) . and (|30|) for the last inequality. 
On the other hand, we have 



We can 
C 



j3 2 h 2 



(#Pr ) 9pr ) i 



(M 1 y, y) = (g Pr , g PT ) + g{M l {gx = (g P r > 9W> + glB A:D R A AR^B^ D g x 



h 2 
1 

7^ 



(g P r ! ffpr 



/I (/) (J) 



□ 



We also need the following lemma. 

Lemma 9 For any w — (w/, pi, wa, wn) £ Wo, i5cw G i?G- 

Proof: We know for any (f 7 , f A , f n ) e Wj ©W A © W n , .g defined by (33 
is in For any u> = (w/, pj, wa, wn) S Woi from the definition of A in 
(fT3|) , there always exists (fj, f A , fn) G W/ © Wa © Wn, such that 



Aw 



f/ 



fA 

fn 



i.e., w — A 1 



f/ 



fA 

fn 



Taking such (f/, fA, fn), 5 defined in (fT2f is Bc»- □ 

The following lemma is an immediate result of Lemmas [8] and |H 

Lemma 10 T/ie space Rq is the same as the range of Be applied on Wq. 

The condition number bound of the preconditioned operator M~ l G is given 
in the following theorem. 

Theorem 4 For all x = (pr, A) G Rm- 1 g> 

C0 1 (Mx,x) < (Gx,x) < ®(H/h) (Mx,x) , 

where §(H/h) is as defined in Lemma\^ f3 as in (0). 



DOMAIN DECOMPOSITION FOR INCOMPRESSIBLE STOKES 



18 



Proof: 

(Gx,x) = x T B c A- 1 B^x = x T B c A- 1 AA- 1 B^x = (l^B^A^B^x 



Since A 1 B^x G Wq and (-, ■) j defines an inner product on VFo, we have 

' u, _ It? \2 

(33) {Gx,x}= max ^ , , L± = max ^^' X ( . 

vGW ,v^O \ v i v )a veW ,v^O /Av,v) 

Lower bound: From Lemma [5J we know that for any given y = (g Pr ,g\) G 
Rg, there exits w e Wq, such that Bqw — y and (^Aw,w^ < jj? (M~ 1 y,y). 
From (|3"3")l . we have 



Since y is arbitrary, using Lemma [6j we have 

(Gx, x) > CP 2 max , ^'^ ; = C/3 2 (Mx, x) . 
yeRa,y^o {M^y^y) 

Upper bound: From ([33]) , Lemmas [3 1 1 01 and |6j we have 

(Gx,x) < $>(H/h) max f lf ■ 
weWo.^o (-M l B c v,B G v) 

= <f(H/h) max =*(H/h)(Mx,x). □ 

Remark 5 From Theorem^and Remark^ we can see that the condition num- 
ber bound of the preconditioned operator M _1 G is independent of the number of 
subdomains when H/h is fixed. If only subdomain corner velocities are chosen 
as coarse level primal variables in the algorithm, the upper eigenvalue bound of 
the preconditioned operator depends on H/h in terms of (H/h) (1 + log (H/h)); if 
both subdomain corner and edge-average velocity degrees of freedom are chosen 
as primal variables, the upper eigenvalue bound grows as H/h. 

Remark 6 With only minor modifications, the algorithm proposed in this paper 
and its analysis apply equally well to the discontinuous pressure case. In that 
situation, p^ and the blocks related to it in ^ can simply be replaced by the 
vector containing subdomain constant pressures and its corresponding blocks, 
respectively. The formulation of the algorithm then follows the same way as 
presented in Section [^J and the same condition number bounds as in Theorem^ 
will be obtained. Numerical experiments of our algorithm for the discontinuous 
pressure case will also be reported in the next section. 
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Remark 7 The same condition number bound has been proved by Kim and 
Lee \14\ \12l with Park] for their FETI-DP algorithms for solving incompressible 
Stokes equations. In their algorithms, discontinuous pressure is considered and 
their approaches do not apply to the continuous pressure case. 

Remark 8 We also note that, no additional coarse level degrees of freedom, 
except those necessary for solving positive definite elliptic problems, are re- 
quired in our algorithm to achieve a scalable convergence rate. For example, 
for two-dimensional problems, it is sufficient to include only the subdomain cor- 
ner velocity degrees of freedom in the coarse level problem. This represents a 
progress compared with earlier work, e.g., [20, 22], where additional continuity 
constraints enforcing the divergence- free conditions on subdomain boundaries 
are required in the coarse level problem. Reduction in the coarse level problem 
size has also been achieved for algorithms discussed in J21 [21 liffl \H\ \12\/ , even 
though discontinuous pressures are considered there. 



7 Numerical experiments 

We consider solving the incompressible Stokes problem ([1]) in the square domain 
SI = [0,1] x [0,1]. Zero Dirichlet boundary condition is used. The right-hand 
side function f is chosen such that the exact solution is 



u 



sin 3 (7rx) sin 2 (7ry) cos(7ry) 
— sin 2 (7ra;) sin 3 (7ry) cos(7rcc) 



and p = x 2 — y 2 



The modified Taylor- Hood mixed finite clement, as shown in Figure [IJ is 
used for the finite element solution. The preconditioned system is solved 
by the CG iteration; the iteration is stopped when the L 2 — norm of the residual 
is reduced by a factor of 10~ 6 . 

Table Q] shows the minimum and maximum eigenvalues of the iteration ma- 
trix M _1 G, and the iteration counts. The coarse level variable space in this 
experiment is spanned by the subdomain corner velocities. We can see from 
Table [1] that the minimum eigenvalue is independent of the mesh size. The 
maximum eigenvalue is independent of the number of subdomains for fixed 
H/h; for fixed number of subdomains, it depends on H/h, presumably in the 
order of (H/h)(l + log (H/h)) as predicted in RemarkO 

For the experiment reported in Table [21 the coarse level variable space is 
spanned by both the subdomain corner velocities and the subdomain edge- 
average velocity components. Even though the edge-average velocity compo- 
nents are not necessary for the analysis, including them in the coarse level 
problem improves the convergence rate, for which the maximum eigenvalue in 
Table [2] grows in the order of H/h, as discussed in Remark [5] 

Tables [3] and [H show the performance of our algorithm for solving the same 
problem, but using a mixed finite element with discontinuous pressure. We 
use a uniform mesh of triangles, shown on the left in Figure [2j the velocity 
finite element space contains the piecewise linear functions on the mesh and the 
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Table 1: Solving (1221) . with only subdomain corner velocities in coarse space. 
H/h (fixed) #sub X min \ max iteration 



8 



4x4 


0.35 


8.92 


21 


8x8 


0.35 


10.07 


28 


16 x 16 


0.35 


10.23 


29 


24 x 24 


0.35 


10.30 


29 


32 x 32 


0.35 


10.33 


29 



#sub (fixed) H/h X min X max iteration 



8x8 



4 


0.30 


4.22 


21 


8 


0.35 


10.07 


28 


16 


0.35 


24.22 


36 


24 


0.35 


40.12 


43 


32 


0.35 


57.15 


50 



pressure is a constant on each union of four triangles as shown on the right in 
the figure. The same mixed finite element has also been used in |22j . 




Figure 2: The mesh and the mixed finite element. 



Comparing Tables Q] and [2] with Tables [3] and [4] we can see that the conver- 
gence rates of our algorithm, using either continuous or discontinuous pressure, 
are quite similar. 
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Table 2: Solving (|2"2"|) . with both subdomain corner and edge-average velocities 
in coarse space. 

H/h (fixed) ^=sub A m j„ \nax iteration 



8 



4x4 


0.36 


4.29 


17 


8x8 


0.36 


5.29 


21 


16 x 16 


0.36 


5.56 


21 


24 x 24 


0.36 


5.61 


21 


32 x 32 


0.36 


5.64 


21 



#sub (fixed) H/h X m in Xmax iteration 



8x8 



4 


0.33 


4.00 


18 


8 


0.36 


5.29 


21 


16 


0.36 


11.63 


26 


24 


0.36 


18.67 


31 


32 


0.36 


26.12 


36 
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